#include "Dwarf.Unmanaged.Hybrid.h"

Solver :: Solver()
{
    length = 0;

    MPI_Comm_rank(MPI_COMM_WORLD, &commRank);           // Retrieve rank of the process
    MPI_Comm_size(MPI_COMM_WORLD, &commSize);           // Retrieve count of processes

    isRootThread = commRank == NUMBER_ROOT_PROCESS;
	W = NULL;
	Wlength = -1;
}

Solver :: ~Solver()
{   
    if(isRootThread)
    {
        for(int i = 0;i < length; i++){
            delete[]complexArray[i];
        }
        delete []complexArray;
    }
	delete W;
}

//creation of work matrix
void Solver :: init()
{
    if(isRootThread)
    {
        this->complexArray = new complex<double>*[length];
        for(int i = 0; i < length; i++)
        {
            this->complexArray[i] = new complex<double>[length];
        }
    }
}

//Transposition of matrix
void Solver :: transpose(complex<double>** matrix)
{
    complex<double> temp;
    for (int i = 0; i < length; i++)
    {
        for (int j = i + 1; j < length; j++)
        {
            temp = matrix[i][j];
            matrix[i][j] = matrix[j][i];
            matrix[j][i] = temp;
        }
    }
}

//1-dimension Fast Fourier Algorithm
void Solver::FFT(complex<double> *complexLine) 
{
    //variables for the FFT
    int bitwidth = (int)(log((double)length)/log(2.0));
    complex<double> w,t,u,v;
    //the complex array splits to real array and complex array 

	if (W == NULL || Wlength != length) 
	{
		Wlength = length;
		W = new complex<double>[length];

		//compute phase weights and reorder for butterflies
		for (int i = 0; i < length; i++)
		{
			W[i].real(cos(-((2* i * M_PI) / length)));
			W[i].imag(sin(-((2* i * M_PI) / length)));
		}
	}

    //reorder for butterflies
    int firstBit = length >> 1;
    int index = 0;
    int reverseIndex = 0;
    for (; ; ) {
        if (reverseIndex < index) {
            complex<double> w = complexLine[index];
            complexLine[index] = complexLine[reverseIndex];
            complexLine[reverseIndex] = w;
        }
        index++;
        if (index == length)
            // Break here, because the code below doesn't terminate for all 1's.
            break;

        // Add 1 to reverse index, going from left to right:
        int carry = firstBit;
        // If adding 1 changed the bit to 0, we have a carry:
        while (((reverseIndex ^= carry) & carry) == 0)
            carry >>= 1;
    }

	//combine and untangle applying phase weights
    //using successive butterflies of 2,4,...,bitwidth 
    for (int i = 1; i <= bitwidth; i++)
    {
        int halfStep = (1 << (i - 1));
        for (int j = 0; j < halfStep; j++)
        {
            int tempShift = (j << (bitwidth - i));
            w = W[tempShift];
            for (int k = j; k < length; k += halfStep * 2)
            {
                u = complexLine[k + halfStep];
                v = complexLine[k];
                t = w * u;
                complexLine[k + halfStep] = v - t;
                complexLine[k] = v + t;
            }
        }
    } 
}

//For the demonstration purposes of the project,
//we have chosen two-dimensional iterative FFT algorithm
void Solver::solve() 
{
    //transform rows
    EvaluateRowsUsingMPI(complexArray);

    /*Transposition to switch rows and columsn*/
    if (isRootThread) transpose(complexArray);

    //transform rows (columns prior to trnspose)
    EvaluateRowsUsingMPI(complexArray);

    /*Transposition to restore rows and columns*/
    if (isRootThread) transpose(complexArray);
}

//Evaluate an FFT on every row in matrix
void Solver::EvaluateRowsUsingMPI(complex<double>** matrix)
{
    MPI_Status st;
    //the temporary complex array
    complex<double> ** tempComplexArray = new complex<double>*[lengthForProc];

    /*send data to slaves*/
    if(this->isRootThread) 
    {   
        int sendTo;
        for(int i = 0; i < lengthForProc; i++)
        {
            tempComplexArray[i] = matrix[i];
        }
        for(int i = lengthForProc; i < length; i++)
        {
            sendTo =  i / lengthForProc;
            if(sendTo > (commSize - 1) )
            {
                sendTo = commSize - 1;
            }
            MPI_Send(matrix[i], length, MPI_DOUBLE_COMPLEX, sendTo, 0, MPI_COMM_WORLD);
        }
    } 
    else
    {   
        for(int i = 0;i < lengthForProc; i++)
        {
            tempComplexArray[i] = new complex<double>[length];
        }
        for(int i = 0; i < lengthForProc; i++)
        {
            MPI_Recv(tempComplexArray[i], length, MPI_DOUBLE_COMPLEX, NUMBER_ROOT_PROCESS, 0, MPI_COMM_WORLD, &st);
        }
    } 
    //end send

    /*Transform the columns(after transposition) */
    #pragma omp parallel for
    for (int j = 0; j < lengthForProc; j++) 
    {
        FFT(tempComplexArray[j]);
    }
    //end transform 
    MPI_Barrier(MPI_COMM_WORLD);

    /*collect results on master */
    if(!this->isRootThread) 
    {
        for(int i = 0;i < lengthForProc; i++)
        {
            MPI_Send(tempComplexArray[i], length, MPI_DOUBLE_COMPLEX, NUMBER_ROOT_PROCESS, 0, MPI_COMM_WORLD);
        }
    } 
    else
    {   
        int receiveFrom;
        for(int i = 0; i < lengthForProc; i++)
        {
            matrix[i] = tempComplexArray[i];
        }
        //collecting of distributed results on master processor
        for(int i = lengthForProc; i < length; i++)
        {
            receiveFrom =  i / lengthForProc;
            if(receiveFrom > (commSize -1) )
            {
                receiveFrom = commSize -1;
            }
            MPI_Recv(matrix[i], length, MPI_DOUBLE_COMPLEX, receiveFrom, 0, MPI_COMM_WORLD,&st);

        }
    }
    //end collection
    if(!this->isRootThread) 
    {
        for(int i = 0; i < lengthForProc; i++)
        {
                delete[]tempComplexArray[i];
        }
    }
    delete []tempComplexArray;
}


